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Abstract 

The current methods to determine the primary energy of ultra-high energy cosmic rays (UHECRs) are different when dealing 
with hadron or photon primaries. The current experiments combine two different techniques, an array of surface detectors and 
fluorescence telescopes. The latter allow an almost calorimetric measurement of the primary energy. Thus, hadron-initiated showers 
detected by both type of detectors are used to calibrate the energy estimator from the surface array (usually the interpolated signal 
at a certain distance from the shower core S (ro)) with the primary energy. On the other hand, this calibration is not feasible when 
searching for photon primaries since no high energy photon has been unambiguously detected so far. Therefore, pure Monte Carlo 
parametrizations are used instead. 

In this work, we present a new method to determine the primary energy of hadron-induced showers in a hybrid experiment based 
on a technique previously developed for photon primaries. It consists on a set of calibration curves that relate the surface energy 
estimator, S (ro), and the depth of maximum development of the shower, X^ax^ obtained from the fluorescence telescopes. Then, 
the primary energy can be determined from pure surface information since S (ro) and the zenith angle of the incoming shower are 
only needed. Considering a mixed sample of ultra-high energy proton and iron primaries and taking into account the reconstruction 
uncertainties and shower to shower fluctuations, we demonstrate that the primary energy may be determined with a systematic 
uncertainty below 1% and resolution around 16% in the energy range from to eV. Several array geometries, the shape 
of the energy error distributions and the uncertainties due to the unknown composition of the primary flux have been analyzed as 
well. 
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1. Introduction 

The energy spectrum of cosmic rays extends by more than 10 
orders of magnitude from below 1 GeV to more than 10^^ eV. 
The energy spectrum follows a power law as E~^, where y is 
around 3.0 in the whole energy range. It is so steep that direct 
measurements are not feasible above 100 TeV. At higher ener¬ 
gies, the properties of the primary cosmic ray are determined 
indirectly from the measurement of the extensive air shower 
(EAS) it produces after colliding with molecules of the atmo¬ 
sphere. 

The highest energy EASs have been traditionally studied us¬ 
ing two different techniques. The first one is based on tele¬ 
scopes that collect the fluorescence light emitted by atmo¬ 
spheric Nitrogen molecules excited by secondary particles of 
the EAS (e.g., Ely’s Eye, HiRes). This allows to determine 
the longitudinal profile of the shower and it is considered to 
be close to a calorimetric measurement of the UHECR primary 
energy. However, fluorescence light can only be observed dur¬ 
ing moonless nights and, consequently, this technique can only 
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be applied to ~ 13% of the incoming events ill]. The sec¬ 
ond technique involves an array of detectors located at ground 
level, mainly scintillators (e.g.. Volcano Ranch, AGASA, KAS- 
CADE) or water Cherenkov tanks (e.g., Haverah Park), whose 
duty cycle is close to 100%. Thus, the lateral distribution of 
secondary particles at ground level can be inferred from the dis¬ 
crete sampling of the shower front. The lateral distribution is 
fitted assuming an appropriate parametrization (called the lat¬ 
eral distribution function, LDE). The interpolated signal at a 
certain optimum distance, S (ro), is used as the energy estima¬ 
tor, which can be related to the primary energy thorough, for 
instance, Monte Carlo (MC) parametrizations. The optimum 
distance, ro, is traditionally fixed for each detector since it is as¬ 
sumed to be only dependent on the array spacing and geometry 
||3] , although some studies suggest the convenience of calculat¬ 
ing the optimum distance for each individual shower taking into 
account its primary energy and direction 11] . 

The current experiments, on the other hand, use hybrid tech¬ 
niques for S (ro) calibration. The Pierre Auger Observatory iH] , 
taking data since 2004, pioneers the simultaneous use of wa¬ 
ter Cherenkov detectors and fluorescence telescopes, while the 
Telescope Array Observatory 11], operating since 2008, com¬ 
bines scintillators and telescopes. Events detected simultane¬ 
ously by both the surface and fluorescence detectors are called 
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hybrid. Hybrid events allow the calibration of S (ro) with pri¬ 
mary cosmic ray energy |@]. Thus, the energy of each event de¬ 
tected by the surface detector alone can be determined almost 
independently of MC simulations. Systematic errors in energy 
estimate are greatly reduced in this way 13, S]. These calibra¬ 
tions assume that the primaries are nuclei and, therefore, they 
cannot be directly applied to photon-initiated showers. In ad¬ 
dition, no photon event has been unambiguously identified up 
to now by any experiment so a proper calibration for photons 
is not possible with this technique. Therefore, each experiment 
relies on MC simulations to infer the primary energy of photon 
events SEimiillil. 

The method used for photon searches by Auger in Ref. iQ 
was first proposed in Ref. (I3l- This method takes into account 
the well-known universality of the electromagnetic component 
of EAs i[il[il[il and the small muon fraction of the photon- 
initiated showers. The calibration curve, that is obtained from 
MC simulations, relates S (ro), the zenith angle of the incoming 
shower, 6, and X^ax- Thus, the primary energy of photon pri- 
maries can be determined with resolution of ~20-25% Eui. 


In this work, we show how to modify that method to be appli¬ 
cable to hadron-initiated showers where the muon component is 
significant, especially, in case of water Cherenkov arrays which 
enhanced their contribution to the total measured signal. The 
additional advantage is that the same method could be used to 
infer the primary energy for both, photon and hadron showers. 
Moreover, in case of hadron-initiated showers the method can 
be calibrated with hybrid events reducing the systematic uncer¬ 
tainties coming from the high energy hadronic models used for 
shower simulations. 


2. Shower and detector simulations 


The simulation of the atmospheric showers is performed with 
the AIRES Monte Carlo program (version 2.8.4a) Q using 
QGSJET-II-03 fl^ as the hadronic interaction model. The in¬ 
put primary energy goes from log(E/eV) = 18.5 to 19.6 in 0.1 
steps. Approximately 2000 events have been simulated per en¬ 
ergy bin for both, proton and iron primaries. The zenith angle 
has been selected following a sine-cosine distribution from 0 
to 60 degrees, while the azimuth angle is uniformly distributed 
from 0 to 360 degrees. X^ax is obtained from these simulations. 

Given the energy, the zenith and azimuth angles of the 
shower, the detector response is simulated with our own code, 
previously tested in Refs. Eollowing the original 

proposal in Ref. iQ, we select a triangular array of cherenkov 
detectors separated 1.5 km and S(ro = 1000m) = 5(1000) as 
the energy estimator. The real core is randomly located inside 
an elementary cell while the reconstructed core position is de¬ 
termined by fiuctuating the real one with a Gaussian function 
whose standard deviation depends on the primary energy, com¬ 
position and the distance between detectors (see Ref. 13] for 
more details). 

The signal collected at each station for a given shower is set 


assuming a true lateral distribution function of the form, 

5(r) = 5(1000) X ^ X ^ , d) 

VoI \^o + g / 

where = 700 m, ro=1000 m, the distance to the shower axis 
r is in meters, 5 (1000) is in VEM (vertical equivalent muons, 
unit for the energy deposited by a vertical muon in a water tank 
13 ]) and)g(^, 5 (1000)) is given by (based on work by T. Schmidt 
et al. ll^ as presented in ll^ ). 

( ab(sQC 6 - 1) ifsec^<1.55 

y5(6>,5(1000)) = i a-^b(sQc6-l) (2) 

[ -r/(sec^ - 1.55)^ ifsec^>1.55 

where a = 2.26-r0.195 log(^), Z? = -0.98, c = 0.37-0.51 sec^-r 
0.30 sec^ 6, d = 1.27 - 0.27 sec^ -r 0.08 sec^ 6, e = c S (1000)^ 
and / = -0.29. 

A realistic 5 (1000) to be used in Eqs. o and (O is obtained 
from. 


E = AiSisf, 

5(1000)(6>) = S3sx[i+Cx-Dx^], (3) 

where v = cos^(6) - cos^(38^). A, B, C and D are constants 
given in Ref. 11231] for QGSJetII-03, iron and proton primaries. 
In addition, shower to shower fiuctuations for each primary are 
emulated by fiuctuating the value from Eq. (O with a Gaussian 
distribution whose standard deviation is taken from Eig. 3 in 
Ref. Q. 

Einally, the signal assigned to each station is fiuctuated using 
a Poissonian distribution whose mean is given by the true LDE. 
We adopt 5= 3.0 VEM and Ssat =1221 VEM as trigger and 
saturation thresholds respectively ||3] • 

Next, the lateral distribution of particles is fitted using a func¬ 
tional form given by. 


log5(r) = ai a 2 



^+G y 

ro + G ) 


(4) 


where the slope of the LDE and the normalization constant are 
free parameters while the core position is fixed in the recon¬ 
structed one. The values of x^jndf are good if at least 3 stations 
are included in the fit, a minimum condition for shower recon¬ 
struction that is fulfilled for almost every event above the energy 
threshold of the detector. Einally, the reconstructed 5 (1000) is 
determined as the interpolated value from the fit at 1000 meters 
from the shower axis. In this method, event by event fiuctu¬ 
ations and reconstruction uncertainties are properly taken into 
account. 

The problem of saturation is common to all surface arrays, 
specially when dealing with high energy vertical showers. The 
consequent lack of detectors close to the core produces large 
uncertainties in the LDE fit and affects the reconstructed 5 (ro). 
The Auger Collaboration, for example, has developed sophis¬ 
ticated algorithms to estimate the signal of a saturated detector 
ll^ . Nevertheless, the analysis of such uncertainties and how 
to minimize them is beyond the scope of the present work so 
saturated events are discarded here. 
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The simulation set has been divided into two samples. In 
each sample, an equal number of proton and iron primaries have 
been mixed for each energy bin. The first sample represents the 
hybrid events and it is used to determine the calibration curves 
as it will be explained in the next Sections. Typical values for 
their reconstruction uncertainties are considered, so their real 
energy, zenith angle and X^ax are fluctuated with Gaussian dis¬ 
tributions whose standard deviations are 15% |@,[^, F 
and 20 g/cm^ iS respectively. The second sample, which rep¬ 
resents data from the surface detector alone, is used for recon¬ 
struction and only their reconstructed S(IOOO) and zenith angle 
are needed. Thus, to simulate the reconstructed 6 we have pro¬ 
ceeded as previously for the calibration sample. Note that the 
method proposed here is obviously also applicable to pure sur¬ 
face arrays but the calibration should be performed with MC 
simulations in this case as, in fact, always occurs in these type 
of experiments. 


3. Method 


The basic idea is that the dependence of S (ro) with energy 
and zenith angle can be factorized as S (ro) = f (0), where a 
is slightly less than 1 (for example 0.95 in Refs. ||l,[l3]). f(6) 
takes into account the longitudinal evolution of the shower so 
it should be a decreasing function of 6 and it depends on the 
slant depth of the shower X = Xgroundl cos(^), where Xground is 
the atmospheric depth at ground level. f{6), as a function of Z, 
behaves similarly to the global profile of the shower. Thus, as 
first approximation, f{6) is a function ofX-Xynax with a similar 
shape to the Gaisser-Hillas function commonly used to describe 
the lon gitu dinal profile. An empirical parametrization is given 
in Ref. lil by. 


5(1000) 

E 


Y AX-100 



(5) 


where AZ = Z - X^ax^ Pi and p 2 are in g/cm^, S (1000) is in 
VEM, the energy is in EeV and po is in VEM/EeV. In the case 
of photon showers, whose muon contribution to the signal could 
be neglected and considering the universality of the electromag¬ 
netic component of EAS, this function results nearly indepen¬ 
dent of the primary energy. In fact, in Ref. Il^ a universal 
parametrization is found with po = lA VEM/EeV, pi = 1000 
g/cm^ and p 2 = 340 g/cm^. 

However, the profile does depend on the primary energy for 
hadron-initiated showers mainly due to the existence of a non- 
negligible muonic component in the showers, which is itself 
a function of the primary energy. Moreover, in case of wa¬ 
ter Cherenkov detectors the sensitivity to muons is enhanced. 
Therefore, despite the fact that the shape of the calibration curve 
is dominated by the photon component of the signal, the muon 
component breaks the universality. In fact, if a unique func¬ 
tion were applied to determine the primary energy of hadron- 
induced showers, it would result in an energy dependent bias, as 
it will be shown later in Sec. lO Therefore, the parameters po, 
Pi and p 2 should be allowed to change with the primary energy 


to correctly reproduce the profile of hadron showers. Then, pi 
tends to be larger than 1000 g/cm^, usually fluctuating around 
3000 g/cm^. In fact, the function is very slightly modified if pi 
is larger than 1000 g/cm^ (the numerator is very close to unity), 
so we have decided to fix pi = 3000 g/cm^. Then, po is the 
maximum of the function and decreases as energy increases. 
Einally, p 2 is related to the width of the function, decreasing 
smoothly as energy increases. 

On the other hand, X^ax in Eq. © can be obtained from its 
average dependence on energy, 

Xjnax = qo-l-qiX log(E/EeV). (6) 

where qo and qi are in g/cm^. 

Therefore, we propose in this work to obtain from the hybrid 
events the next calibration curves: 

(A) the global curve using all these events following Eq. ©, 

(B) a set of curves, one for each energy bin, following Eq. 

in order to account properly for the muonic component of 
the showers. It is important that these curves do not cross 
and that the statistics is good enough to assure a good fit 
near their maximum, 

(C) the Xmax evolution as a function of energy following 

Eq.®, 

where the parameters needed are S (1000), the zenith angle, the 
reconstructed energy and Xynax^ They could be obtained from 
the standard fluorescence reconstruction and the LDE fit. These 
curves, obtained from simulations, are shown in Eig. [T] (more 
details in Sec. SHi. 

Then, given a pure surface event, its energy could be deter¬ 
mined from the reconstructed S (1000) and the zenith angle of 
the incoming shower. Both can be obtained from the LDE fit 
and the geometrical reconstruction respectively. The procedure 
is as follows: 

(1) Using an initial estimation of the primary energy (5 ,10 or 
30 EeV), Xynax is obtained with (C). As it will be shown 
later, the reconstructed energy do not depend on this choice. 

(2) Xmax is used to get an energy estimation using (A). 

(3) This energy is used to get X^ax again with (C). 

(4) Steps (2) and (3) are repeated until the difference be¬ 
tween two consecutive energies converges to a stable value 
(AEjE < 10“^). Around 3-5 iterations are required. 

At this step, a reconstructed energy is obtained, but as explained 
before, the non-universality of (A) for hadron primaries intro¬ 
duce a significant bias that must be corrected. 

(5) The previous reconstructed energy is used to select the 
nearest calibration curve from the set (B). 

(6) Steps (1) to (4) are repeated using the new curve from (B) 
instead of (A), so a new reconstructed energy is obtained. 

(7) If the nearest curve from (B) to the new reconstructed en¬ 
ergy is the same as before, the process finishes. Otherwise, 
(5) and (6) are repeated. Only 2 or 3 different calibration 
curves are usually needed. 
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Following this procedure the energy bias is corrected as it will 
be shown later in Sec. lO 

In a real experiment, the events with a large error in the re¬ 
constructed zenith angle should be analyzed carefully or even 
rejected since the process could not converge. In fact, they are 
mostly saturated events with energy very close to the threshold 
of the detector, or events with cores very close to the border 
of the array or to a detector that is not working, so they are in 
general already rejected by the quality cuts usually imposed for 
data analysis. 


4. Results 

4.1. Calibration 

Independent sets of calibration curves can be, and indeed 
were, obtained for each primary. Nevertheless, we only show 
here the results for an equal mixture of proton and iron pri¬ 
maries. The global curve and the fits for each energy bin 
(Alog(£’/eV) = 0.1) are presented in Fig. [T]a and[T]b respec¬ 
tively. The higher the energy, the lower the curve. Parameters 
po and p 2 are free while p\ is fixed at pi = 3000 g/cm^ as pre¬ 
viously explained, po and p 2 smoothly decrease as energy in¬ 
creases. The evolution of X^ax as a function of energy is shown 
in Fig.[T]c, where the medians of the distributions are fitted tak¬ 
ing into account the corresponding confidence levels shown in 
the figure. Note that the first and the last energy bins are not 
included since both are unavoidably affected by the limited en¬ 
ergy range of the simulation set. Therefore, despite the fact 
that their corresponding curves are shown in Fig.[T]b, neither of 
them is shown in the remaining plots. 

4.2. Energy reconstruction 

We use the reconstruction sample, which is statistically inde¬ 
pendent from the calibration, in order to test the method. Given 
the reconstructed 5(1000) and zenith angle of each event, its 
energy is determined. As an example. Fig. [2] shows the itera¬ 
tive process for a typical event. Solid line represents the global 
calibration curve while the dashed line is the one used in the 
last step of the process. It can be seen that the latter is much 
closer to the real position of the event (red star), improving 
significantly the energy reconstruction. It was verified that the 
convergence point is independent of the path followed, so the 
reconstructed energy does not depend on the initial value used 
to start the iteration as previously mentioned. 

Fig.E shows the error in the reconstructed energy. The two 
bins at the edges have been rejected for the same reasons as 
before. Fig. [3la shows that the energy dependent bias resul¬ 
tant from the use of only the global calibration curve is greatly 
reduced by employing a set of calibration curves for discrete 
energies, although a residual error <1% still remains. Fig. Ob 
also shows the error for proton and iron primaries. The error 
is almost energy independent for both, each primary taken sep¬ 
arately and also for the mixed composition sample. Note also 
that the error bars, which can be considered as the resolution 
of the method, are also slightly energy dependent. The energy 



Xground/cos0 - Xmax [g/cm 2] 




logio(E/EeV) 


Figure 1: Calibration curves for a mixed sample of iron and proton primaries. 
The global fit (a) and a different fit for each energy bin (b) are shown. In (b) 
the curve is lower as energy increases. The evolution of X^ax as a function of 
energy is shown in (c). See Sec[3for details. 
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Figure 2: An example of the reconstruction method: first, the global calibra¬ 
tion curve (solid line) allows to get a first estimation of the primary energy by 
iterating between Eq. e and itself (process not shown). This energy is used 
to select a new calibration curve (dashed line) repeating the iterative process 
with the new curve instead of the global curve. The path followed during the 
iterative process is shown in the Figure. As can be seen, the convergence point 
is the same independently of the initial value used to start the iteration (5, 10, 
30 EeV) and it is very close to the real position of the event (red star). In this 
example, the error in the final reconstructed energy is -5% while it would be 
larger than -i-30% if only the global curve were used. 


error distribution for both primaries taken together is approx¬ 
imately Gaussian and the resolution of the method is around 
16%, as shown in Fig. |4l In the case of Iron and proton taken 
separately, the distributions are also nearly Gaussian and the 
resolution is around 13.5% and 14.5% respectively. 

It is important to note that the method is quite sensitive to the 
^max VS. log(£’/EeV) calibration curve. In fact, if the parameter 
qo in Eq. ® is modified by +10%, an energy error of +7% will 
be produced. This effect is not so important for qi, since if it 
changes ±10%, a negligible error of +1% will be obtained. 

4.3. Gaussianity of the energy error distributions 

Arguably, it is desirable that the errors in energy reconstruc¬ 
tion follow a Gaussian distribution. Gaussian errors, for exam¬ 
ple, are easier to handle and understand when applying decon¬ 
volution techniques for determination of the spectrum while as¬ 
suring that there are no asymmetries or long tails, which further 
reduces the danger of the limited energy range of the detector 
and biases associated with a rapidly changing spectral index. 
An example of these undesirable effects can be seen in Ref. ^ . 

In case of a Gaussian distribution, the ratio between the high 
and low parts of the 68% and 95% confidence levels (C.L.) 
should be 1 since the distribution is symmetric. In addition, the 
ratio between the 95% and 68% C.L., for both the high and low 
parts, should be 2 since they represent the values for 2cr and Icr 
respectively. Eig.[5]shows both ratios for the energy error distri¬ 
bution obtained in each energy bin. It has been calculated using 
the bootstrap technique that consists on resampling the original 
distribution a large number of times and calculating these ratios 
for each new sample. The median an 68% C.L. for each ratio 
are shown in the plots. It can be seen that the error distribu¬ 
tions are nearly Gaussian with a very good symmetry and the 
absence of long tails. 
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Figure 3: Energy error as a function of primary energy. The points and the error 
bars are the median and the region of 68% probability, respectively, (a): Using 
only a global calibration (red squares) or using a full set of calibration curves 
(black circles) for a mixture of iron and proton primaries, (b): Proton (blue 
triangles) and iron (red squares) taken separately and together (black circles) 
are shown. 



Figure 4: Energy error distribution for the whole sample where equal number 
of iron and proton primaries in each energy bin are selected. 
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Figure 5: Analysis of the symmetry (top) and the tails (bottom) of the energy 
error distributions. See text for details. 




Figure 6: Energy error as a function of the proton fraction of the reconstructed 
sample using three different sets of calibration curves: the first one obtained 
with proton showers, the second one from only iron events and, finally, taking 
both type of primaries together (called mixed here) as shown in Sec. l4Jl The 
points and the error bars are the median and the confidence levels at 68 and 
95% respectively. 


4.4. Energy bias as a function of the primary composition 

We have determined the energy error as a function of the 
composition of the reconstruction sample. Since the energy 
error and resolution are almost energy independent for every 
primary (Fig. O, we have mixed all the events in this analysis 
regardless of their energy. We have used 100 samples with 100 
events each. Proton and Iron primaries are randomly selected 
such the proton fraction varies from 0 to 1 in 0.1 steps. The en¬ 
ergy error as a function of the proton fraction is shown in Fig.0 
The errors are also shown for the case that the calibration curves 
would have been obtained for each primary separately (we call 
them proton and iron calibrations in the figure). As expected, 
the error is negligible if a sample with proton fraction of 0.5 is 
reconstructed with the calibration from Fig. [T] since equal frac¬ 
tion of proton and iron events were used to get this calibration 
curves. Note also that if each primary is reconstructed with 
its own calibration the error is also very close to zero. In the 
unrealistic scenario that an extremely pure composition sample 
were reconstructed with the calibration obtained from a mixed 
composition fiux, the absolute value of the incurred error would 
be < 10%. 

4.5. Different array sizes and geometries 

The robustness of the method has been tested by varying the 
geometry of the array, i.e. the shape of the unitary cell and the 
distance between detectors. 

First, an array with detectors arranged in a squared layout, 
with a separation of 1200 m is analyzed. Such a geometry is 
characteristic, for example, of the Telescope Array (TA) obser¬ 
vatory. The applicability of the method to a scintillator array as 
TA is not analyzed here since it would require a different study 
as a consequence of the different response of scintillators to the 
muonic part of the shower compared to Cherenkov tanks. As 
mentioned, this Section focuses on its application to Cherenkov 
tank arrays with different geometries. The calibration curves 
and the error in the reconstructed energy are shown in Fig.|7]for 
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Figure 7: Calibration curves (top) and error in the reconstructed energy (bot¬ 
tom) for the Telescope Array geometry. The points and the error bars are the 
median and the confidence levels at 68%. 


a mixed sample. Note that S (800) has been selected as the en¬ 
ergy estimator 0 . The error is energy independent and smaller 
than 1%, while the resolution is -16%. 

Second, a smaller triangular array is considered. A geome¬ 
try representative of the Auger Infill is selected, i.e. a triangular 
grid with 750 m spacing between detectors. The Infill is es¬ 
sentially designed to explore the region below the ankle and to 
reach full efficiency above -10^^-^ eV ll^ . The fraction of sat¬ 
urated events is very large above 10^^ eV, so we only analyzed 
here the interval from 10^^-^ to 10^^ eV. In the energy range of 
the experiment, the average optimum distance is rp = 450 m, so 
S (450) is selected as the energy estimator ll^ . Then, the error 
in the reconstructed energy is around -5% as shown in Fig.[8l 
The error comes from the lack of events close to the maximum 
of the calibration curves which, as mentioned previously, is a 
key point in order to get a reliable fit and calibration. Since the 
optimum distance increases with primary energy 13], the aver¬ 
age ro is higher for the energy range analyzed here. Thus, the 
error could be minimized if a larger rp were selected. For ex¬ 
ample, using S (750) the error would be negligible (Fig. [8]). 


5. Conclusions 

An iterative method, previously developed to infer the pri¬ 
mary energy of photon-induced showers in pure surface arrays, 
has been modified to be applicable to hadron-initiated showers 
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Figure 8: Error in the reconstructed energy in case of a triangular array whose 
distance between detectors is 750 m. Two different energy estimators are used, 
*S(450) and *S(750). Points and error bars are the median and 68% confidence 
levels. See text for details. 


and tested assuming a hybrid experiment. The method is based 
on a set of calibration curves that relate the surface energy es¬ 
timator, S (ro), and the depth of maximum development of the 
shower, thanks to hybrid events. In pure surface arrays, a 
similar procedure could be performed but the calibration should 
rely on Monte Carlo (MC) parametrizations. However, it is im¬ 
portant to note that MC parametrizations could be affected by 
the fact that the simulations do not reproduce properly the avail¬ 
able experimental data S, a major advantage of the hybrid 
experiments. 

The original method is based on the well-known universality 
of the electromagnetic component of the showers and the small 
number of muons produced in the development of photon cas¬ 
cades. However, the significant fraction of the muon compo¬ 
nent for hadron-initiated showers, breaks the universality and 
makes necessary to implement a full set of calibration curves 
depending on the primary energy. 

Our own simulation program of the detector response has 
been used. Shower to shower fluctuations and reconstruction 
uncertainties (core position, LDF fit and signal fluctuations) 
have been implemented. Primary energy and zenith angles 
go from to eV and from 0 to 60 degrees respec¬ 
tively. Several array geometries varying the shape of the uni¬ 
tary cell (triangular and square) and the distance between de¬ 
tectors (1500, 1200 and 750 m) have been studied. Considering 
a mixed sample of proton and iron primaries, the energy is de¬ 
termined with a negligible error and resolution around 16% in 
the full energy range analyzed. 

Obviously, the same composition is expected for hybrid 
(used for the calibration) and pure surface events, so the en¬ 
ergy of the latter could be determined with error close to zero. 
However, in the extreme scenario that the reconstructed events 
present a pure composition, the energy error could achieve 
-10%. This could be considered as the maximum uncertainty 
of the method due to the unknown composition of the primary 
flux. 

The energy error distributions are nearly Gaussian, an impor¬ 
tant point to get a more reliable reconstruction of the shape and 
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position of rapidly varying spectral features since they are eas¬ 
ier to manage when applying deconvolution techniques in the 
spectrum determination. 
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